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INTRODUCTION 


Picatinny  Arsenal  Technical  Report  #4831  entitled  "Dynamic  Thermal 
Modeling  of  the  Behavior  of  a Nuclear  Warhead  Section  Adaption  Kit  in 
an  Accidental  Fire  Environment"  treated  the  problem  of,  a safety  analysis 
based  on  dynamic  thermal  modeling.  The  purpose  of  this  report  is  to 
provide  potential  users  with  a description  of  the  thermal  modeling 
method  itself  and  the  details  of  the  computer  programming  method 
developed. 

It  was  pointed  out  in  Picatinny  Arsenal  Technical  Report  #4831 
that  unknown  effective  emissivities  and  convective  film  coefficients 
would  be  determined  empirically  by  performing  one  or  a few  simulated 
fire  tests  run  under  closely  controlled  conditions, , and  then  fitting 
the  model  temperatures  versus  time  to  the  test  thermocouple  results. 

The  general  lumped-constant  (lumped-parameter)  method  for  solving 
complex  dynamic  engineering  problems  is  well  known  (Reference  1). 

The  method  of  analysis  developed  here  for  the  thermal  behavior  of  a 
warhead  section  adaption  kit  is  of  this  type. 

MODELING  METHOD 

The  specific  device  which  was  modeled  is  the  Spartan  Warhead 
Section  Adaption  Kit.  It  is  comprised  of  some  26  different  types  of 
mechanical  parts  and  electrical  components.  It  is  three-dimensional, 
involving  several  different  materials,  complex  geometric  shapes,  and 
multipath  heat  flow  patterns.  Time  is  an  essential  part  of  the  problem 
so  the  model  must  be  a dynamic  one. 

The  first  step  in  the  lumped-constant  method  is  to  divide  the 
device  to  be  modeled  into  a number  of  pieces  called  elements.  The 
temperature  of  each  element  is  the  temperature  of  a single  point 
within  the  element.  In  the  model,  there  are  no  temperatures  other 
than  that  of  each  of  the  elements.  Heat  is  transferred  from  one 
element  to  another  by  conduction,  radiation  and/or  convection.  The 
model  is  a network  made  of  the  elemental  points  with  heat  paths 
between  some  pairs  of  points.  The  temperature  of  an  element  rises  at 
a rate  which  depends  on  the  net  heat  flow  to  that  element  and  its 
heat  capacity.  All  of  the  heat  capacity  associated  with  an  element 
is  assumed  to  be  concentrated  at  the  elemental  point. 

The  simulation  modeled  here  is  based  on  a mental  image  of  a possible 
accident  scenario.  If,  before  attaching  the  first  stage  to  a Spartan 
missile,  an  organic  fuel  fire  were  to  start,  the  Warhead  Section 
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Adaption  Kit  cover  could  be  exposed  to  the  flames  but  its  base, 
being  mounted  against  the  warhead,  would  not.  Under  worst  case  con- 
ditions such  a fire  is  considered  equivalent  to  exposure  of  the  cover 
to  an  185CP  F black  body  source  of  heat.  Average  organic  fire  con- 
ditions would  correspond  to  about  145(f  F. 

In  this  situation  the  cover  will  quickly  heat  up  and  radiate  heat 
to  components  inside.  There  will  also  be  conduction  to  the  base  and 
through  the  upright  brackets  to  the  deck.  Convection  will  also  convey 
heat  from  the  cover  to  cooler  surfaces  inside.  Dense  smoke  will  be 
generated  when  plastic  insulation  and  paint  overheats  and  this  may 
affect  heat  flow.  It  is  expected  that  some  of  the  aluminum  and  mag- 
nesium parts  may  melt  since  their  melting  points  are  well  below  the 
heat  source  temperature.  If  the  exposed  part  of  the  deck  melts,  a 
new  radiation  path  is  created  through  the  hole  to  the  base  below. 

Each  of  the  effects  described  above  was  accounted  for  in  the  model. 

The  Spartan  Warhead  Section  Adaption  Kit  was  divided  into  74 
elements.  (The  number  of  elements  chosen  is  arbitrary,  but  a large 
number  becomes  unwieldy  whereas  a small  number  may  not  include  enough 
detail  to  accurately  distinguish  temperatures  between  the  crucial 
points.  The  number  of  elements  is  limited  by  the  number  of  represen- 
tation statements  required,  which  cannot  exceed  685.  The  program 
described  here  has  about  600.)  The  way  in  which  the  division  was 
done  was  dictated  by  the  results  desired.  In  the  case  treated  here, 
the  temperature  versus  time  was  required  for  eight  critical  points 
in  order  to  study  several  "thermal  races". 

In  a region  where  two  such  points  are  close  together,  the  division 
has  been  made  fairly  fine  while  other  regions  are  coarsely  sectioned. 

The  propellant  activated  battery  has  two  critical  points,  a gas 
generator  which  turns  it  on,  and  a thermal  release  connector  which 
disconnects  it  electrically  from  the  rest  of  the  circuit.  (Actually, 
there  are  two  gas  generators  and  two  thermal  release  connectors,  but 
attention  is  concentrated  on  the  set  which  is  known  to  be  most  critical.) 
Because  of  the  importance  of  the  thermal  race  between  these  two  points 
and  the  geometry  in  this  region,  the  battery  has  been  divided  into 
27  elements. 

Elements  are  identified  by  a two  digit  number.  Elements  of  the 
base  plate  are  numbered  01  to  45,  of  the  deck  plate  50  to  64,  uprights 
(which  support  the  deck  plate  on  the  base  plate)  80  to  89.  The 
warhead  thermal  release  connector  90  and  the  cover  97  to  99.  The 
shroud  which  was  used  as  a heat  source  surrounding  the  whole  adaption 
kit  is  called  element  SH. 
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In  the  computer  program  to  be  described  TC22  is  a variable  name, 
the  value  of  which  at  any  given  time  is  the  temperature  of  thermocouple 
number  22.  Similarly,  T56  is  a variable  name  the  value  of  which  is  the 
temperature  of  element  number  56.  All  other  thermocouples  and  elements 
have  corresponding  variable  names. 

Figures  1,  2 and  3 show  in  outline  the  important  parts  of  the 
Spartan  Warhead  Section  Adaption  Kit  which  have  been  modeled.  Figure 
1 is  a side  view  of  a section  through  the  center  of  the  adaption  kit. 

The  cover  and  its  three  elements  can  be  seen.  Element  99  is  the  part 
of  the  cover  below  the  flange.  Element  98  is  that  part  of  the  cover 
above  the  flange  to  the  level  of  the  deck  plate.  Element  97  is  the 
rest  of  the  cover.  Elements  98  and  99,  instead  of  representing  a 
point  within  the  cover,  are  represented  by  a line  going  around  the 
cylinder  halfway  between  the  top  and  bottom  of  that  part  of  the  cover. 
Figure  1 also  shows  the  base  plate  at  the  bottom  and  the  deck  plate 
part  way  between  it  and  the  top  of  the  cover. 

Figure  2 is  a plan  view  of  the  base  plate  with  its  major  components. 
The  base  plate  is  divided  into  12  elements  numbered  as  shown  (01  through 
12). 

The  battery  (element  numbers  19-45)  occupies  the  major  part  of  the 
lower  half  of  the  figure  and  a Burst  Delay  Timer  (element  numbers 
13-18)  the  upper  part.  The  prefixes  to  the  element  numbers  (used 
only  on  the  diagram)  designate: 


T - 

Top 

B - 

Bottom 

U - 

Upper 

L - 

Lower 

M - 

Match  (Propellant  Element) 

S - 

Side 

The  elements  80  through  39  are  the  upright  supports  for  the  deck 
plate.  The  outer  supports  (80-86)  are  fastened  to  the  cover  with 
screws.  Element  90  is  the  warhead  connector. 

Figure  3 is  a plan  view  of  the  deck  plate.  Mounted  on  it  are  the 
IC  Box  and  two  arm-safe  devices  (ASD's).  Element  56  is  ASD  Nr.  1 and 
corresponds  to  thermocouple  TC22  as  shown  in  Figure  3.  The  correspondence 
between  other  elements  and  thermocouple  numbers  in  the  heat  test  is 
given  in  the  table  on  page  17. 
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MATHEMATICAL  APPROACH 


1.  General . 

The  dynamic  simulation  implemented  here  is  represented  by  74 
simultaneous  differential  equations: 


^ = 10'^  QnFn/Vn  “F/second  (1) 


in  which: 

Tn  = temperature  of  n'th  element  in  *^F 
t = time  in  seconds 

On  = heat  flux  toward  n'th  ela.ient  i.i  e^^Vsec 
Vn  = volume  of  n'th  element  in  cubic  inches 

n n - 

Cn  = specific  heat  of  n'th  element  in  Cal/g/  C 
Dn  = density  of  n'th  element  in  g/cc 

Each  Qn  is  the  sum  of  several  terms  which  represent  the  heat  flux 
toward  element  n from  each  other  element.  These  heat  flux  terms  are 
of  three  types; 

conduction  QC  « fj  (T^  - T^) 
radiation  QR  = (Tj^  - T^^) 

convection  (similar  to  conduction) 

The  conduction  terms  involve  elements  in  contact.  There  is  one 
such  term  for  each  pair  of  elements  in  contact. 

Radiation  terms  involve  elements  between  which  there  is  a direct 
path  for  radiation.  The  only  significant  ones  Involve  one  very  hot 
element  and  appreciable  areas  facing  one  another.  Since  these  criteria 
also  apply  to  convection,  the  same  element  pairs  are  used  for  both. 

2.  Terminology. 

The  rate  of  heat  flow  (in  erg/sec)  is  designated  Q In  this  report. 
This  heat  flux  from  element  01  to  element  02,  when  It  results  from 
conduction,  is  called  QC0102.  For  radiation  from  99  to  01,  the  flux 
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is  QR9901.  The  total  flux  toward  an  element  01  is  called  QOl.  It  is 
the  algebraic  sum  of  all  terms  which  involve  that  element.  For  example 

QOl  = QR9901  - QC0102  (Z) 

The  second  term  is  negative  because  the  order  of  the  element 
numbers  indicates  heat  flow  away  from  01.  (Eq.  (2),  which  illustrates 
the  terminology,  is  incomplete.  The  actual  expression  includes  12  more 
conduction  terms  and  two  more  radiation  terms.  The  complete  equation 
is  given  in  the  Appendix.) 

There  is  one  such  equation  for  each  element  of  the  model. 

3.  Conduction  - (Ref  2). 

Elements  in  contact  have  heat  conduction  terms  between  them  (See 
Figure  4).  The  heat  flux  expression  is  of  the  form: 

QC0102  = KC*(T01-T02)/R  (See  Note  below)  (3) 

where  R is  the  thermal  resistance  and  KC  is  a constant  for  unit 
conversion  (1.411  here).  The  resistance  R in  (3)  is  the  sum  of  three 
terms,  that  within  element  01  to  the  boundary,  the  contact  resistance 
and  that  within  02  from  the  boundary  to  the  elemental  point.  This 
can  be  expressed  as: 

R = D01B/(K01*A01)  + D02B/(K02*A02)  + R0102C  (4) 

DOIB  = Distance  heat  must  travel  within  element  01  to  the  boundary 
within  element  01  in  inches 

AOl  = Cross-sectional  area  for  heat  flow  in  element  01  (Sq  inches) 

KOI  = Thermal  conductivity  of  element  01  (erg/cm-sec-deg  C) 

D02B,  A02,  K02  similarly 
R0102C  = Contact  resistance 


NOTE:  Equations  in  this  text  and  in  the  computer  program  are  written 

in  the  format  of  Fortran  in  which: 

* signifies  multiplication 
/ signifies  division 
**  signifies  exponentiation 

parentheses  (multi-nested,  if  necessary)  are  used  to  resolve 
ambiguities. 

3.48E-5  stands  for  3.48  x 10'^ 


In  cases  where  the  two  elements  are  part  of  the  same  material,  or 
welded  together,  the  contact  resistance  is  assumed  to  be  zero.  The 
distance  within  element  01  is  the  distance  from  the  defined  point  to 
the  boundary  with  02.  The  cross-sectional  area  for  conduction  is  an 
estimated  average  area  between  the  point  and  the  boundary.  For  an 
element  which  is  a rectangular  parallelepiped,  it  is  the  full  cross- 
sectional  normal  to  the  direction  of  heat  flow.  This  is  justified  on 
the  basis  that  the  conductive  heat  flow  is  not  toward  the  point  but 
rather  is  practically  one  dimensional.  The  temperature  gradient  is 
nearly  uniform  for  the  whole  region.  Since  all  the  heat  conduction 
paths  (except  the  three  to  the  battery  explosive  actuator)  are  metallic, 
it  becomes  evident  when  analyzing  the  actual  test  data  that  all  of  the 
gradients  are  quite  small  so  the  linear  approximation  is  believed  to 
be  reasonable. 

If  we  substitute  (4)  into  (3)  and  use  width  times  thickness  for 
cross-sectional  area,  we  get: 

QC0102  = 1.411*(T01-T02)/(D01B/(K01*W01*TH01)  + 002B/(K02*W02*TH02)  + 
R0102C)  (5) 

which  is  another  of  the  basic  equations  used  in  the  computer  program. 

There  is  one  such  equation  for  each  pair  of  elements  in  contact. 

4.  Radiation  - (Ref  3 and  4). 

In  the  terms  we  are  using  in  this  analysis,  radiative  heat  flow 
is  of  the  form: 

QR9901  = KR  *((T99  + 459.7)**4  - (TOl  + 459.7)**4)  * G * AOl  * V9901 

(6) 

KR  = Stefan-Boltzmann  constant  and  a factor  for  units  (3.48E-5  Here) 
AOl  = Area  of  element  01  in  square  inches 
V9901  = Effective  view  factor 

G = 1/(1/E01+1/E99-1)  (See  Ref  4,  P74,  Eq.  (4.7))  (7) 

EOl  and  E99  are  effective  emissivities  and  both  elements  are 
assumed  to  be  gray  bodies. 

G has  a value  between  zero  and  one.  It  is  like  a composite  emissivity. 
(Ref  3,  P76,  Eq  at  bottom  of  page  is  similar.  In  this  case  one  body 
is  a total  enclosure  so  its  emissivity  is  1 in  which  case  (7)  reduces 
to  G = E). 
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V9901  is  the  view  factor  or  configuration  factor.  It  has  a value 
between  0 and  1 which  represents  the  fraction  of  element  99  which 
is  "seen"  by  01  and  vice  versa. 

Equation  (7)  assumes  the  absorptivity  is  equal  to  the  emissivity 
and  all  reflection  is  diffuse.  It  takes  into  account  radiation, 
absorption,  reflection  and  reradiation. 

Because  of  the  fourth  power  of  absolute  temperatures,  it  is 
evident  that  the  magnitude  of  radiative  heat  flux  is  very  dependent 
on  temperature.  Only  terms  of  which  one  of  the  elements  is  at  high 
temperature  need  be  considered.  The  only  ones  considered  here  are 
those  in  which  one  of  the  cover  elements  is  included.  By  the  time 
any  other  element  reaches  a comparable  temperature,  some  melting  occurs. 
The  most  exposed  areas  inside  the  cover  are  the  deck  plate  and  the 
IC  Box,  which  melt  at  1110  and  1150  degrees,  respectively. 

5.  Convection. 

The  convective  heat  flux  is  expressed  as: 

Q = hA  At  (Ref  5,  page  3)  (8) 

where  h = coefficient  of  heat  transfer 

A = Area  of  Suv'face 
At  = Temperature  Difference 

The  whole  subject  of  convective  heat  transfer  is  very  complex  and 
quantitative  results  are  difficult  to  compute.  But  for  the  present 
purpose,  it  has  been  assumed  that  convective  heat  flux  is  proportional 
to  the  temperature  difference  with  the  constant  of  proportionality 
determined  empirically.  It  is  also  assumed  that  only  those  element 
pairs  involved  in  radiative  heat  flow  may  also  include  the  effects 
of  convection.  This  simplifies  the  coding  considerably,  without,  it 
is  believed,  detracting  from  the  utility  of  the  program. 

SYSTEM  SIMULATOR 


1.  General . 

The  digital  computer  program  used  as  a system  simulator  was  written 
in  CSSL  3 (Continuous  System  Simulator  Language).  A detailed  under- 
standing of  the  program  requires  familiarity  with  CSSL  3,  which  is 
described  in  the  User's  Manual,  (Ref  6),  and  with  FORTRAN.  A brief 
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description  of  the  CSSL  3 features  which  have  been  used  will  be  given 
here. 

Statements  may  be  either  in  FORTRAN  or  CSSL  3 syntax.  Most 
statements  may  be  placed  in  any  order  within  the  program.  A sort 
routine  arranges  the  order  of  replacement  statements  so  that  all 
variables  are  defined  before  they  are  used.  The  programmer  can  force 
a fixed  order  for  a group  of  statements  by  the  use  of  a PROCEDURAL. 

All  statements  (CSSL  3 or  FORTRAN)  are  written  in  free  format  in 
which  any  column,  from  1 through  72,  may  be  used.  A statement  may  be 
continued  onto  the  next  line  by  three  periods  before  the  end  of  the 
line.  Blanks  may  be  added  freely. 

At  the  beginning  of  a statement,  a variable  name  (alphameric) 
followed  by  one  period  denotes  a logical  control  variable.  If  that 
variable  is  true  at  execution  time,  the  statement  will  be  executed; 
otherwise  not.  A statement  which  starts  with  the  word  CONSTANT 
indicates  that  the  variable  names  which  follow  have  the  value  given 
which  is  invariable  with  time  (the  only  independent  variable  in  the 
program)  but  which  may  be  given  a different  value  at  execution  time 
for  a different  run. 

2.  Tabular  Functions. 

A function  may  be  inputted  to  the  program  as  a series  of  points 
called  a TABLE.  For  example,  TABLE  SHROUD,  1,  14,  etc.,  means  that 
a Table  function  has  the  name  SHROUD,  one  dependent  variable  and  14 
pairs  of  values  (14  points)  in  which  the  next  14  values  are  of  the 
independent  variable  in  order  of  increasing  value  and  the  remaining 
14  values  are  of  the  corresponding  dependent  variable.  Interpolations 
between  points  are  linear.  Extrapolations  are  linear  from  the  last 
two  points.  This  function  is  used  in  the  statement: 

NX.TSH  = SHROUD  (T) 

which  indicates  that,  for  the  present  value  of  T,  TSH  is  to  be  set 
to  the  corresponding  value  of  the  SHROUD  function  (if  NX  is  TRUE). 

Data  for  this  table  was  obtained  during  the  fire  simulation  test 
by  digitally  recording  thermocouple  voltages  on  computer  tape.  Tabular 
function  SHROUD  was  used  as  the  driving  function  of  temperature  versus 
time  for  the  differential  equation  set.  Similarly,  the  eight  other  sets 
of  thermocouple  data  were  inputted  as  table  functions  to  facilitate 
comparisons  with  model  results. 


description  of  the  CSSL  3 features  which  have  been  used  will  be  given 
here. 

Statements  may  be  either  in  FORTRAN  or  CSSL  3 syntax.  Most 
statements  may  be  placed  in  any  order  within  the  program.  A sort 
routine  arranges  the  order  of  replacement  statements  so  that  all 
variables  are  defined  before  they  are  used.  The  programmer  can  force 
a fixed  order  for  a group  of  statements  by  the  use  of  a PROCEDURAL. 

All  statements  (CSSL  3 or  FORTRAN)  are  written  in  free  format  in 
which  any  column,  from  1 through  72,  may  be  used.  A statement  may  be 
continued  onto  the  next  line  by  three  periods  before  the  end  of  the 
line.  Blanks  may  be  added  freely. 

At  the  beginning  of  a statement,  a variable  name  (alphameric) 
followed  by  one  period  denotes  a logical  control  variable.  If  that 
variable  is  true  at  execution  time,  the  statement  will  be  executed; 
otherwise  not.  A statement  which  starts  with  the  word  CONSTANT 
indicates  that  the  variable  names  which  follow  have  the  value  given 
which  is  invariable  with  time  (the  only  independent  variable  in  the 
program)  but  which  may  be  given  a different  value  at  execution  time 
for  a different  run. 

2.  Tabular  Functions. 

A function  may  be  inputted  to  the  program  as  a series  of  points 
called  a TABLE.  For  example,  TABLE  SHROUD,  1,  14,  etc.,  means  that 
a Table  function  has  the  name  SHROUD,  one  dependent  variable  and  14 
pairs  of  values  (14  points)  in  which  the  next  14  values  are  of  the 
independent  variable  in  order  of  increasing  value  and  the  remaining 
14  values  are  of  the  corresponding  dependent  variable.  Interpolations 
between  points  are  linear.  Extrapolations  are  linear  from  the  last 
two  points.  This  function  is  used  in  the  statement: 

NX.TSH  * SHROUD  (T) 

which  indicates  that,  for  the  present  value  of  T,  TSH  is  to  be  set 
to  the  corresponding  value  of  the  SHROUD  function  (if  NX  is  TRUE). 

Data  for  this  table  was  obtained  during  the  fire  simulation  test 
by  digitally  recording  thermocouple  voltages  on  computer  tape.  Tabular 
function  SHROUD  was  used  as  the  driving  function  of  temperature  versus 
time  for  the  differential  equation  set.  Similarly,  the  eight  other  sets 
of  thermocouple  data  were  inputted  as  table  functions  to  facilitate 
comparisons  with  model  results. 


3.  Macros. 

CSSL  3 allows  the  usage  of  a MACRO  function.  A MACRO  is  created 
by  a statement  MACRO  MACRO  followed  by  a MACRO  name  and  a dummy 
variable  name.  It  ends  with  a statement  MACRO  END.  In  between,  the 
MACRO  is  defined  with  the  dummy  variable  subscripted  with  digits 
indicating  the  order  of  the  parameters  in  the  call.  An  example  from 
our  program  will  illustrate. 

MACRO  MACRO  COND  Q 

Q(l)  = 1.411E+7  * (Q(2)-Q(3))/(Q(4)  etc 

MACRO  END 


if  called  by  - 

COND  QC0102,  TOl,  T02,  4.9  etc. 
will  result  in  - 

QC102  = 1.411E+7*(T01-T02)/(4.9,-*-  etc 
being  inserted  into  the  program. 

Three  CSSl  3 MACRO'S  were  written.  The  first  called  COND  and 
used  in  the  example  above  produces  heat  flux  values  at  each  time 
step  due  to  conduction  between  elements  pairs  using  the  form  of  Eq.  (5). 
The  second,  called,  RAD,  similarly  produces  heat  flux  due  to  radiation 
and  convection  using  the  form  of  Eqs.  (6)  and  (8).  The  third,  called 
OT,  produces  equations  of  the  form  of  Eq.  (1). 

4.  Procedurals. 

A PROCEDURAL  BLOCK  is  used  to  force  a fixed  order  of  execution 
for  some  statements.  This  begins  with,  for  example,  PROCEDURAL  (M52, 
M58,  M19,  M56,)  etc.,  and  ends  with  END.  This  block  is  sorted  so 
that  all  variables  to  the  right  of  the  equal  sign  come  before  the 
block  and  those  to  the  left  come  afterwards.  This  particular 
PROCEDURAL  simulates  the  melting  of  elements  52,  58,  19  and  56  (see 
the  complete  program  listing  in  Appendix).  For  these  four  elements, 
in  addition  to  the  variables  T52,  etc.,  two  more  were  set  up  for  each 
element. 


0 = no 

1 = yes 


Element  Temp 

Element  Melting  Point 

Melted 

T52 

MP52 

M52 

T58 

MP58 

MSB 

T19 

MP19 

M19 

T56 

MP56 

M56 

Each  M- variable  such  as  M52  is  set  to  zero  unless  the  corresponding 
temperature  T52  is  equal  to  or  greater  than  the  melting  point  of 
element  b2,  MP52  in  which  case  M52  is  set  to  one.  M52  is  used  as  a 
multiplier  for  terms  which  are  to  be  added  when  melting  has  occurred, 
such  as  radiation  from  element  97  to  element  01  through  the  void. 

The  expression  (1.  - M52)  is  used  as  a multiplier  for  the  heat  flux 
to  element  52  and  also  for  each  other  term  which  should  be  dropped  at 
this  time  such  as  conduction  to  adjacent  elements. 

Another  PROCEDURAL  BLOCK  is  used  to  switch  the  SHROUD  temperature 
versus  time  from  the  tabular  function  to  an  exponential  function  with 
controllable  parameters. 

The  third  PROCEDURAL  BLOCK  is  for  the  simulation  of  the  effect 
of  smoke  on  radiative  heat  flux.  (The  reason  for  introducing  this 
effect  is  given  in  the  subsequent  section  "Fitting  the  Model  to  Heat 
Test  Results".) 

A fourth  PROCEDURAL  allows  switching  off  the  heat  input  to 
elements  98  and/or  99  to  simulate  the  situation  when  the  whole  cover 
is  not  exposed  to  heat. 

5.  Data  Inputted. 

Data  is  inputted  to  the  program  in  the  form  of  constants,  MACRO 
call  parameters  or  table  values. 

Data  inputted  as  program  constants  include: 

Coefficient  of  Convective  Heat  Transfer 

View  Factor 

Emissivity 

Melting  Points 

Contact  Resistances 

Specific  Heat  i Density  via  Factor  F (See  Mathematical  Approach 

General ) 
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Data  inputted  as  MACRO  call  parameters  are: 

MACRO  PARAfiETERS 

COND  Distance  for  conductive  heat  flow 

Width 

Thickness  For  cross-sectional  area  for  conduction 
Conductivity 
RAD  Emissivity 

Area  Exposed  to  Radiation 
View  Factor 

DT  Specific  Heat  & Density  via  Factor  F 

Volume  of  Element 

The  only  data  inputted  by  table  values  are  temperatures  versus 
time  obtained  from  physical  testing.  Nine  such  tables  were  used  in 
the  program. 

6.  Fitting  the  Model  to  Heat  Test  Results. 

Once  all  of  the  known  physical  and  thermal  data  have  been  inputted 
to  the  model,  the  next  step  is  to  find  values  for  the  various  unknowns 
(such  as  convective  heat  coefficients)  which  make  the  model  correspond 
to  the  test  results  reasonably  well.  There  are  eight  thermocouple 
temperature  versus  time  curves  derived  from  the  test,  each  of  which 
correspond  to  an  element  of  the  model.  Values  are  arbitrarily  chosen 
for  each  unknown  constant  and  the  eight  curves  produced  by  the  model 
are  compared  with  the  test  curves.  On  the  first  guess  Tiost  of  the 
model  curves  differ  considerably  from  the  test  curves.  Because  some 
model  curves  are  influenced  more  by  some  constants  than  others,  it  is 
not  difficult  to  make  the  general  trend  of  the  curves  match  the  test 
curves  in  about  20  runs.  (Incidentally,  these  runs  do  not  require  the 
complete  CSSL  cycle.  The  first  part  produces  and  compiles  a FORTRAN 
program  which  takes  substantial  computer  time.  But  this  does  not 
have  to  be  redone  if  the  program  itself  is  not  changed.  It  is  only 
necessary  to  change  some  of  the  “constant"  values  and  this  can  be 
done  at  run  time.  Each  run  itself  then  takes  very  little  time.) 

After  getting  the  curves  to  generally  fit  overall  temperature  rise, 
the  next  step  is  to  attempt  to  get  the  shape  of  the  curves  right  by 
balancing  the  contributions  due  to  radiation,  convection  and  conduction. 
All  curves  were  matched  reasonably  well  in  this  way  except  for  thermo- 
couple TC22  in  the  test  which  corresponds  to  element  T56  in  the  model. 
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This  thermocouple  is  on  the  deck  plate  not  covered  by  the  IC  Box  and 
is  thus  exposed  to  radiation  and  convection  from  the  two  upper  parts 
of  the  cover.  The  conduction  from  the  cover  through  the  deck  plate 
itself  allows  the  use  of  only  one  arbitrary  unknown  constant  which 
is  the  cover  to  deck  plate  contact  resistance.  But,  even  with  this 
at  zero,  the  contribution  to  temperature  rise  from  conduction  is  very 
small.  This  means  that  the  heating  of  the  deck  plate  is  mainly  by 
radiation  and/or  convection.  To  find  out  how  much  of  each  may  be 
required,  two  sets  of  runs  were  made.  In  the  first  set,  convection 
was  set  at  zero  and  radiation  increased  and  decreased  until  the  latest 
significant  point  (the  occurrence  of  melting)  is  matched.  The  resultant 
curve  is  shown  on  Figure  5.  In  the  second  set  radiation  was  eliminated 
to  the  deck  plate  and  with  only  convection  the  end  point  was  again 
matched.  That  curve  is  also  shown  on  Figure  5,  together  with  the 
curve  for  the  test  result.  The  problem  is  evident  in  the  Figure. 
Combinations  of  radiation  and  convection  can  only  give  curves  between 
the  two  extremes  whereas  the  test  curve  is  outside.  Therefore,  the 
actual  phenomenon  is  not  merely  heating  by  radiation  and  convection 
(with  conduction,  too,  of  course).  Something  else  is  going  on. 

During  the  test  it  was  observed  that  a large  volume  of  dense 
smoke  was  evolved.  It  was  felt  that  this  smoke  may  affect  heat  flow. 

If  so,  it  would  have  an  effect  which  changes  with  time.  This  could 
affect  the  heating  curve  in  the  way  observed,  with  the  heating  rate 
higher  in  the  beginning  and  decreasing  with  time.  To  simulate  this 
a smoke  function  (S)  was  introduced.  The  precise  manner  in  which 
smoke  is  generated  is,  of  course,  unknown.  It  seemed  reasonable  to 
assume  that  smoke  will  begin  to  form  when  the  temperature  is  high 
enough  for  some  materials  to  start  burning.  Four  hundred  degrees 
Fahrenheit  was  chosen  as  the  temperature  at  which  smoking  begins. 

The  value  chosen  is  probably  not  critical  since  the  temperatures 
rise  rapidly.  It  was  also  felt  that  the  rate  of  smoke  generation 
increases  more  rapidly  than  proportional  to  temperature  so  a square 
function  was  used.  The  amount  of  smoke  present  is  cumulative  so  the 
rate-of-smoke  generation  DS  is  integrated  outside  the  PROCEDURAL  to 
give  S.  The  effective  net  emissivity  GV  is  then  made  a simple  inverse 
function  of  S. 

Using  this  smoke  function  the  arbitrary  constants  Kl,  K2  and  K3 
were  adjusted  to  make  T56  match  TC22  while  keeping  all  other  curves 
reasonably  well-matched.  Figures  6 through  13  are  the  eight  pairs  of 
curves  for  the  model  output  compared  with  the  thermocouple  readings 
when  the  empirical  curve-matching  was  concluded. 
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7.  Computer  Program. 

The  Appendix  has  a complete  source  listing  of  both  the  main  program, 
which  must  be  run  once  each  time  the  model  is  changed,  and  of  a typical 
parameter-alteration  run. 

The  main  program  is  in  three  contiguous  parts.  The  first  part 
(lines  100-230)  is  the  control  deck  for  the  SCOPE  operating  system. 

The  second  part  (lines  250-7030)  is  the  CSSL  3 source  deck.  The 
third  part  (lines  7050-7180)  is  the  CSSL  3 run  time  control  deck. 

Execution  proceeds  in  several  steps.  First,  CSSL  3 reads  the 
CSSL  3 source  deck,  sorts  the  statements,  makes  MACRO  substitutions, 
interprets  statements  which  are  in  CSSL  3 format  and  produces  a FORTRAN 
source  program  on  a file  called  TAPE7.  This  is  then  compiled  with  the 
machine  code  written  on  a file  named  LGO.  Compiled  subroutines  on  a 
permanent  file  are  loaded  into  central  memory  along  with  LGO  for  execution 
as  the  last  step.  Total  time  for  execution  of  the  main  program  is 
typically  250  CPU  seconds. 

A parameter  alteration  run  is  illustrated  by  the  second  source 
listing  in  the  Appendix.  Any  defined  constants  in  the  main  program 
(see  lines  430-520)  may  have  new  values  assigned  for  this  type  of 
run.  The  time  step  and  ending  time  may  be  reset  by  redefining  Cl  and 
TE  respectively.  The  constant  XX,  initially  set  to  zero,  makes  the 
shroud  temperature  equal  to  the  table  function  which  is  the  measured 
shroud  thermocouple  values  as  a function  of  time.  If  reset  to  1 
an  exponential  temperature  rise  is  used  with  a maximum  temperature  and 
time  constant  of  TMAX  and  TAU  respectively.  Similarly,  the  constants 
X98  and  TBC  can  be  used  as  switches  to  turn  on  or  off  the  radiation 
to  elements  T93  and  T99  respectively.  The  values  of  any  variable 
versus  time  may  be  outputted  and/or  plotted  during  a parameter 
alteration  run.  This  allows  detailed  study  of  heat  flow  within  the 
model.  Execution  of  a parameter  alteration  run  uses  5 to  10  seconds 
CPU  time. 


USING  THE  MODEL 


There  are  several  ways  in  which  the  model  can  now  be  used.  One 
of  these  is  to  assume  a different  heat  input.  The  most  interesting 
heat  input  is  a shroud  temperature  which  is  a fast  exponential  rise  to 
a steady  1850°F.  This  corresponds  to  a worst-case  organic  fuel  fire. 

The  result  of  using  such  an  input  on  T52  is  shown  in  Figure  14.  As 
expected,  of  course,  the  rate  of  heating  is  greater  than  was  experienced 
in  the  test,  but  in  addition,  details  such  as  the  time  of  melting  can 
probably  best  be  estimated  by  the  use  of  a model  like  this.  The  change 
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to  the  program  which  is  required  to  get  this  result  is  to  input  a 
value  for  XX  of  1 instead  of  0.  This  makes  the  logical  variable  X 
true  through  PROCEDURAL  (X.NX  = XX)  which  acts  as  a switch  to  change 
the  definition  of  TSH,  the  temperature  of  the  shroud  versus  time. 

With  XX  = 0,  TSH  is  the  SHROUD  tabular  function  of  thermocouple  read- 
ings whereas  when  XX  = 1,  TSH  equals  an  exponential  function  of  time. 

TMAX  is  set  to  1850°and  TAU  is  10  (seconds)  within  the  program. 

Figure  15  was  gotten  in  the  same  way  but  the  maximum  shroud  temper- 
ature (Tf4AX)  was  set  to  1450“F  which  is  considered  to  be  typical  of 
an  organic  fuel  fire  rather  than  worst-case. 

Another  configuration  which  was  considered  is  1850“F  heat  only  to 
the  top  of  the  cover,  which  could  occur  if  the  missile  were  lying  on 
the  ground  with  a fuel  fire  nearby.  The  resulting  T52  is  shown  in 
Figure  16.  To  simulate  this  condition  two  switches  were  reset  in  the 
program.  The  constant  TBC  was  reset  to  2000,  and  X98  was  reset  to  1. 

TBC  is  in  the  program  as  the  time  at  which  the  missile  ballistic  case 
melts.  Prior  to  TBC  there  is  no  heat  input  to  element  T99  because  of 
the  shielding  effect  of  the  ballistic  case.  By  setting  TBC  = 0,  we 
have  heat  input  to  T99  all  the  time  whereas  by  setting  it  high  the 
heat  never  gets  inputted  to  T99.  Constant  X98  is  used  only  to  allow 
switching  off  the  heat  input  to  element  T98. 

Temperature  T52  for  the  case  of  heat  input  to  the  whole  top  and 
with  a maximum  input  temperature  of  550(f  F to  correspond  to  a propellant 
fire  is  given  in  Figure  17.  In  a case  like  this,  which  differs  so 
widely  from  the  conditions  under  which  the  model  parameters  were  derived, 
the  results  must  be  interpreted  with  a great  deal  of  caution.  Never- 
theless, such  results  may  be  better  than  a simplified  calculation  or 
estimate. 

Another  use  for  the  model  is  in  the  analysis  of  heat  flow  paths. 

As  an  example,  consider  how  the  battery  match  gets  heated. 

T45  gets  heated  by  conduction  from  T32,  T34,  and  T42,  which  three 
elements  form  a bottom  corner  of  the  battery  case.  Since  any  variable 
can  be  outputted  (not  only  the  state  variables  used  in  the  integration 
routines)  all  that  is  needed  is  to  output  some  of  the  QC's  and  QR's 
and  compare  these  to  find  the  main  heat  flow  patterns. 

Figure  18  is  a diagram  of  all  major  (>  10^  erg/sec)  heat  flow 
paths  feeding  elements  32.  34  and  42  at  the  time  240  seconds.  It  is 
seen  that  elements  32  and  34  are  receiving  more  heat  by  radiation  than 
conduction.  Element  42  receives  heat  only  by  conduction.  There  are 
several  more  heat  flow  paths,  such  as  radiative  and  conductive  flux 
from  98  to  99  to  01  and  02  which  are  not  shown.  Also,  everything 
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changes  with  time.  All  of  this  information  is  available  to  use,  keeping 
in  mind  that  the  interpretation  must  be  a reasonable  consequence  of  the 
data  put  into  the  model. 

Knowledge  of  the  heat  flow  paths  allows  the  designer  to  incorporate 
conductors  and  insulators  to  obtain  desired  performance.  It  may  be 
desirable  to  insulate  a battery  gas  generator  for  example  while  increasing 
the  conduction  to  the  thermal  release  connector. 
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Table  1 


Designations  of  test  thermocouples  (TC) 
and  computer  model  temperatures  (T) 


Thermocouple  Model 


TC 

25 

T50 

TC 

24 

T52 

TC 

22 

T56 

TC 

27 

T19 

TC 

12 

T45 

TC 

15 

T58 

TC 

21 

T90 

TC 

26 

T07 

Location 

Deck,  under  IC  Box 
Deck,  exposed  part 
Arm  Safe  Device  Number  1 
Battery  TRC 

Battery  explosive  gas  generator 
IC  Box  top 
WHD  TRC 

Base  between  BDT  and  battery 
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Fig  1 Cover 


Fig  3 
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T19  (SOLID)  t TC87  (DOTTED)  OS  TINE 


Fig  9 Computer  model  T19  and  test  thermocouple  TC  27 


T45  (SOLID)  & TC12  (DOTTED)  US  TIHE 
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Fig  10  Computer  model  T45  and  test  thermocouple  TC  12 


T90  (SOLID)  & TC21  (DOTTED)  US  TinE 
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Fig  12  Computer  model  T90  and  test  thermocouple  TC  21 


Te?  (SOLID)  4 TC86  (DOTTED)  US  TINE 
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Fig  13  Computer  model  T07  and  test  thermocouple  TC  26 


T52  (SOLID)  i TC84  (DOTTED)  US  TIME 


Computer  model  heat  source  at  14 SOT 


52  (SOLID)  k TC24  (DOTTED)  US  TIHE 


Computer  model  heat  source  at  1800°F  (to  top  only) 


Tsa  (SOLID)  & TCa4  (DOTTED)  US  TIME 
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Fig  17  Computer  model  heat  source  at  5500°F 


Numbers  in  boxes  are  element  numbers. 
Other  numbers  are  heat  values  in  millions 
of  ergs  per  second  at  240  seconds . 


Fig  18  Heat  paths  to  element  45 
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